Role of the on-site pinning potential in establishing quasi-steady-state conditions of 

heat transport in finite quantum systems 
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We study the transport of energy in a finite linear harmonic chain by solving the Heisenberg 
equation of motion, as well as by using nonequilibrium Green's functions to verify our results. The 
initial state of the system consists of two separate and finite linear chains that are in their respective 
equilibriums at different temperatures. The chains are then abruptly attached to form a composite 
chain. The time evolution of the current from just after switch-on to the transient regime and then 
to later times is determined numerically. We expect the current to approach a steady-state value 
at later times. Surprisingly, this is possible only if a nonzero quadratic on-site pinning potential 
is applied to each particle in the chain. If there is no on-site potential a recurrent phenomenon 
appears when the time scale is longer than the traveling time of sound to make a round trip from 
the midpoint to a chain edge and then back. Analytic expressions for the transient and steady-state 
currents are derived to further elucidate the role of the on-site potential. 

PACS numbers: 05.70.Ln,44.10.+i,63.22.-m 
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I. INTRODUCTION 

In the study of thermal transport in quantum systems, 
the conventional setup examined is a scattering region 
sandwiched between two infinite leads acting as large 
heat reservoirs [l|. The leads and the scattering region 
are initially prepared to be at their respective thermal 
equilibriums in canonical distributions and then are cou- 
pled using an adiabatic switch-on. The long-time steady- 
state current through the scattering region can be calcu- 
lated using the Landauer formula where the transmis- 
sion coefficient can be determined using nonequilibrium 
Green's function techniques lj. This approach, however, 
is not applicable to situations where the system is under- 
going dynamical changes and is far from the steady-state 
regime. An example would be the transient behavior of a 
system where the coupling between the scattering region 
and the leads is not weak and is switched on abruptly ■ 
In addition, the method does not provide any information 
on how the steady state is dynamically approached from 
the initial configuration of the system. In this work we 
examine the time evolution of the thermal current from 
the transient regime just after an abrupt switch-on to the 
intermediate-time quasi-steady-state regime and then to 
the long-time recurrent regime, in a system consisting 
only of two finite leads that are coupled abruptly at time 
t = 0. We examine if, when, and how the onset of the 
steady-state occurs. To determine the current, we solve 
the quantum equations of motion of the linear system. 
In addition, a supplementary second method we employ 
to verify our results is an approach using nonequilibrium 
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Green's functions that takes into account the full time 
evolution of the system Q. 

1*0 1*= 1*0 1*0 ¥*o 1*0 
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(b) combined system 

FIG. 1. (Color online) An illustration of how two finite in- 
dependent systems at temperature TL and Tr, containing Nh 
and particles, are abruptly combined at t — into a com- 
posite system containing A^l + A^r particles. The interparticle 
spring constant is k and the on-site spring constant is fcrj. 

We determine the thermal current in a phonon system 
consisting of two linear chains abruptly attached together 
at time t — 0. Shown in Fig. fflja) are the two sepa- 
rate linear chains. The interparticle spring constant is k 
and the on-site spring constant is fcn. The left and right 
leads contain Nr. and Nr sites and have temperature Tr. 
and Tr, respectively. The chains are initially in their 
respective thermal equilibrium; i.e., they are initially at- 
tached to heat baths so that they acquire their corre- 
sponding temperatures and then the baths are subse- 
quently detached. The chains therefore satisfy canonical 
distributions. They also follow fixed boundary conditions 
wherein particles at the left and right edges are attached 
to fixed walls. At time t = the chains are abruptly cou- 
pled with spring constant k, as shown in Fig. QIb). The 
composite chain also satisfies fixed boundary conditions. 
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We then determine the energy current flowing through 
the newly formed inter leads coupling. Since the chains 
are finite, we cannot use the Landauer formula approach 
to calculate the current. However, results from the Lan- 
dauer formula can be used for comparison to the infinitely 
large and long-time limit results of the approaches we de- 
scribe. 



The probability of extracting work in a system that is 
weakly coupled to finite heat baths is found to follow a 
power law [3|. In this paper, in comparison, we exam- 
ine the energy current flowing between two finite chains 
whose coupling is not weak. We find that the presence of 
a quadratic on-site pinning potential is necessary in es- 
tablishing a steady-state current. Previous studies have 
found that the on-site pinning potential plays an impor- 
tant role in how the steady-state phonon current depends 
on the system size in disordered harmonic systems 
In this work we further investigate the role of the on- 
site potential on the dynamics of the system and find 
that, without a quadratic on-site pinning potential, the 
current exhibits oscillatory behavior that does not disap- 
pear even for long times. Furthermore, we find that the 
time period of oscillation of the current is proportional 
to the sum of the length of the finite chains. Studies in 
systems consisting of classical harmonic oscillators and 
classical particles indicate that the Poincare recurrence 
time, i.e., the time for a specific phase-space configura- 
tion of the system, or a configuration that very closely 
resembles it, to reappear increases exponentially with the 
number of degrees of freedom [j| @ . Recurrence of the 
wave function is also found to occur in quantum systems 
with discrete energy spectrum [7[ and systems that are 
periodically driven Our work extends these studies 
on quantum systems by examining a physical observable, 
i.e., the energy current, in a finite system to see if it 
displays recurrent behavior and to determine conditions 
that dictate the appearance of recurrence. In Sec.lIVIour 
results show that the presence of an on-site potential is 
crucial to determine whether the energy current exhibits 
a recurrent and oscillatory behavior or a behavior that 
decays and settles to a quasi-steady-state value. 



In electron systems, the transport of electrons between 
two finite leads has previously been studied ||. When 
there is a potential bias between the leads, it is found 
that a quasi-steady-state current with a finite lifetime 
appears, even when there are no dissipative effects like 
electron-electron and electron-ion interactions. In this 
paper, we supplement this electron transport study by 
examining the transport of phonons within a finite sys- 
tem. We determine the dynamical behavior of the ther- 
mal current and find that having on-site pinning poten- 
tials is necessary in establishing the quasi-steady-state 
current. 



II. THE EIGENMODE APPROACH 

We model the left and right chains by the harmonic 
Hamiltonian 

H a = ±(p a ) T p a + ±(u a ) T K a u a , a = L, R, (1) 

where u a is a column vector whose elements are the 
rcnormalized displacements of the sites in chain a, p a 
is the conjugate momentum, and K a is the tridiago- 
nal spring constant matrix consisting of 2k + fco along 
the diagonal and — fc along the two off-diagonals, k is 
the nearest-neighbor spring constant and fco is the spring 
constant of the on-site pinning potential. The spring po- 
tentials we consider in this work are all quadratic. The 
equations of motion of particles in the left chain are 

u L = p L and ii L =p L = -K L u L - V hR u R , (2) 

where V LR is the coupling matrix between the left and 
right chains. Particles in the right chain satisfy a similar 
set of equations of motion. The current flowing out of 
the left chain can be calculated from the chain's average 
rate of decrease in energy, 

I L = - (ij L ) = (p L (t) T V LR u R (t)) 

jk 

where the equations of motion in Eq. ([2]) are used. The 
average is taken with respect to the initial product state 
density matrix. To calculate the energy current out of the 
left chain, therefore, we need to determine the dynamics 
of u R (t) andp L (i). 

Consider an isolated finite linear chain consisting of 
./V sites and satisfying fixed boundary conditions. Since 
it is isolated, the equations of motion of its constituent 
particles are 

ii=p and il = p=~Ku, (4) 

where K is the spring constant matrix of the whole chain. 
The solution is 

u(t) — cos(VKt) uq H ^= sm(yfKt) p , 

yK ' (5) 

p(t) = -VK sin(VKt) u Q + cosi^Kt) po, 

where uq and po are determined from the initial condi- 
tion. Since the system is linear, quantum Heisenberg 
operators and classical variables have identical solutions. 
To make the solution numerically amenable, we construct 
a transformation, using eigenmodes, that diagonalizes 
the matrix K. Let ti" be the nth eigenmode of K, i.e., 
K u n — f^u™, where is the eigenvalue associated 
with u n , i.e., 

^ = 2fc(l-cosg n )+fco, (6) 
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where q n = irn/(N + 1) and n — 1, 2, ■ ■ • ,N. Because 
of the fixed boundary conditions, the jth element of the 
eigenmode is u™ = A n sin(q n j), where A n can be fixed 
by an appropriate normalization. Let S be a matrix con- 
sisting of the eigenmodes of K, i.e., S — (u , u 2 , . . . , u w ). 
We can then fix A n by normalizing 5S" T = S T S = 1. We 
get 



(7) 



The matrix K is diagonalized by the similarity transfor- 
mation 

s t ks = di&g{tii,n 2 2 ,...,n 2 N ) = n 2 , (8) 

where the right-hand side is a diagonal matrix consisting 
of elements f2 2 , ft 2 ,, . . ., The coefficients involving 

K in Eq. ([5]) can now be calculated with the aid of the 
above similarity transformation. We are then left with 
the unknowns uq and pa, and their correlations. To de- 
termine these unknowns, we write Uq in normal mode 
coordinates, i.e., uo — SQ, where Q contains the nor- 
mal modes. In a harmonic oscillator with Hamiltonian 
Hho = Ylj Klj^jdj + 1/2), the jth normal mode is 



(9) 



Pj = Qj 



where aj(t) — <Zj-(0)e _ 3 ' and its complex conjugate are 
the time-dependent lowering and raising ladder opera- 
tors. These ladder operators satisfy expectation values 

(a)a^j = S 3k f 3 , ^4) = 5jk(fj + 1), and (a]alj = 

(ajCLk) = 0, where fj — (cxp(MXj /fcgT) — 1) is the 
Bose-Einstein distribution function. Using these expec- 
tation values, we get 



(QjPk) = j s jk , 
(PjQk) = —Sjk, 



(10) 



(p i fl fe ) = n i -(2/ J + i)<y ifc . 

Using S and the above expectation values, we can write 
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For the composite chain in Fig. [ljb) , we label the sites 
consecutively from 1 for the leftmost site to ATl + ATr for 
the rightmost site. Thus, the labels of the sites where 
the newly formed interleads spring appears are ATl and 
+ 1. In Eq. ([3]) therefore, the indices are j = A^ and 
k = Al + 1. The expectation value in Eq. ([3]), using the 
solutions in Eq. ([3]), can then be written as 



,R 
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where 



D j A k (p^ k )+D j B k (p^ k ), 



(12) 



(13) 



and a — L,R. The (up) and (pu) terms cancel exactly 
and only the (uu) and (pp) correlations contribute to the 
currents. Note also that correlation between the left and 
right regions vanishes and each region has its own initial 
temperature T a . The coefficients are the symmetrized 
version of those in Eq. ([5]): 

N 



Ak — ^ Sjv L +i,m cos (0 OT t) S km , 
D a sm(tt m t) 



m— 1 
N 



-5, 



(14) 



^ ' SjVl, m^m sin (Q m t) Sj 

m— 1 
JV 

Dj = Sn^,™ cos (O m f) 5j m , 



where A~ = ATl + A/r. The current can therefore be cal- 
culated using Eqs. ©, (UU), ([13]), (HI]), and flU}. In ad- 
dition, Eq. (|12p can be further simplified into two terms. 
One of the terms is proportional to (/l — /r). This term 
leads to the steady-state value in the limit A" — > oo. The 
other term is proportional to (/l+/r+1/2) and produces 
the transient current. See the Appendix. 

The expression for I L (t) in Eq. (|3]) is the current flow- 
ing through the rightmost site of the left lead, i.e., at the 
site labeled Al- At any site n in the composite lead, a 
general expression for the current flowing through it can 
be written as 



I n (t) = V nm ( Pn {t)u m (t)) 

y ] ^ nm.jk + ^ ] ^ nm,jk 
j,k=l j,fc=iVL + l 



(15) 



where U„ m is the coupling between sites n and m = n + 1 
and 



Tnmjk — CnjA m k (uojUofe) + C n jB m k (uQjPok) 

+ D nj A mk (PojU^k) + D nj B mk (PojPok) 



(16) 



4 



C n j, and D n j are the generalized 



where a = L, R. Equation ([16]) is in the same form as 
Eq. (HU)) but generalized to include the current flowing 
through any site n in the composite chain. Similarly, the 
coefficients A mk , B mk 
forms of those in Eq. (fT4j) . 



III. NONEQUILIBRIUM GREEN'S FUNCTIONS 
APPROACH 
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An alternative approach is to use nonequilibrium 
Green's functions (NEGF) techniques to calculate the 
time-dependent energy current Q . The retarded Green's 
function for a finite collection of harmonic oscillators in 
a chain in equilibrium can be written as 10] 



g r {t) = -S 0(t) 



sin (Tit) T 



(17) 



where S is the matrix in Eq. ([SJ and fl is the diagonal ma- 
trix of the square root of the eigenvalues in Eq. ©. The 
advanced, lesser, and greater equilibrium Green's func- 
tions can then be determined from the above retarded 
Green's function. From these equilibrium Green's func- 
tions, the time-dependent nonequilibrium Green's func- 
tions and the energy current can be calculated following 
the procedure described in Ref. Q. With the eigenmode 
approach in Sec. Inland the NEGF approach discussed in 
this section, we can therefore calculate the energy current 
in two different and independent ways. We find that both 
methods produce the exact same results, up to double- 
precision accuracy. 



IV. NUMERICAL RESULTS AND DISCUSSION 

The energy current can be calculated using either the 
eigenmode approach discussed in Sec.|TI]or the NEGF ap- 
proach described in Sec. lIIII However, the NEGF calcula- 
tion is computationally intensive because of the presence 
of several multiple integrals whose numerical convergence 
must be carefully determined. In contrast, calculations in 
the eigenmode approach involve, as the most complicated 
part, fast and straightforward matrix manipulations. Af- 
ter verifying that our results are the same for several sets 
of data from both approaches, we proceed and acquire 
most of our data using the eigenmode approach. 

We first determine the steady-state heat current and 
conductance of an infinite linear chain using the Lan- 
dauer formula with unit transmission for frequencies 
within the phonon band [ll|. The value of the nearest- 
neighbor spring constant we use is k = 1 eV/A 2 u. We 
examine the steady-state current for on-site pinning po- 
tentials k = 0, 0.01, 0.1, and 1 eV/A 2 u. Shown in Fig. M 
are the plots of the steady-state thermal current and con- 
ductance as a function of the temperature. We use these 
results for comparison to the quasi-steady-state values 
arising in finite leads. 



FIG. 2. (Color online) Plots of (a) the steady-state cur- 
rent and (b) the conductance as functions of the average 
temperature between the two leads and with a tempera- 
ture bias of a = 0.1. The values of the on-site spring con- 
stant are fco = (black line), fco = 0.01 eV/(A 2 u) (red 
dash line), fco = 0.1 eV/(A 2 u) (green dash-dot line), and 
fco = 1 eV/(A 2 u) (blue dash-double dot line). The inter- 
particle spring constant is k = 1 eV/(A 2 u). 
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FIG. 3. (Color online) Plots of the current as a function of 
time when the left and right chains are finite with lengths 
(a) N L = N R = 100 and (b) N L = N R = 50. The average 
temperatures between the leads are T = 10 K (black line), 
T — 100 K [red (gray) line], and T = 300 K [blue (dark gray) 
line]. The temperature bias between the leads is a = 0.1. 
There is no on-site potential, i.e., fco = 0. 




FIG. 4. (Color online) A close up view of the current when 
the average temperature of the leads is T = 10 K. The darker 
(blue) line is when TVl = -/Vr = 100 while the lighter (orange) 
line is when Nl = -/Vr = 50. The dashed line is the value 
of the steady-state current calculated using the Landauer for- 
mula. 



The dynamical energy currents are determined for var- 
ious lengths and initial temperatures of the finite leads 
and the strengths of the on-site spring potential. Shown 
in Fig. |3] are plots of the current flowing out of the 
left lead when there is no on-site spring potential. The 
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length of the leads are — Ar = 100 in Fig. |3fa) and 
N L = N R = 50 in Fig. Mb). The left lead initially has 
temperature Tl = (1 + a)T, where a = 0.1, while the 
right lead has initial temperature Tr = (1 — a)T '. The 
plots in Fig.[3]correspond to T = 10 K, 100 K, and 300 K. 
The leads are attached at time t — 0. Just after the 
switch-on, the current initially shoots down towards neg- 
ative values. Since we are calculating the current flowing 
out of the left lead, a negative current value implies that 
energy is flowing into the left lead. For the right lead, 
we also find the same initial negative current values, i.e., 
energy is also flowing into the right lead. This result 
is the same as that found in systems with semi-infinite 
leads 0j. The reason is because an external input en- 
ergy is required to make the interleads coupling between 
the two previously unattached leads. After the leads are 
connected, the extra energy that is applied to make the 
connection then flows into the leads, thus initially pro- 
ducing negative current values just after the switch-on. 
This phenomenon also occurs for the leads-center system 
and is explained by a small t expansion 



12] 



On longer time scales, the current rises almost linearly 
and then suddenly drops. The overall behavior is roughly 
periodic with a period proportional to the full length of 
the chain. A ringing current can also be observed in 
Fig. M This ringing current also appears in systems con- 
taining infinite leads [2]. Such a behavior has also been 
observed previously in NEGF calculations in electronic 
transport jl3| . In Fig. M notice that the current does 
not stabilize to a quasi-steady-state value, even when the 
leads have the low average temperature T = 10 K, as 
shown in Fig. [4l 




150 200 250 
time [10" 14 s] 

FIG. 5. (Color online) Plots of the current as a function of 
time when the on-site spring constant is ko = 0.001 eV/A 2 u. 
The left and right chains have lengths (a) TVl = Nr_ = 100 
and (b) JVl = Nr = 50. The average temperatures between 
the chains are T = 10 K (black line), T = 100 K [red (gray) 
line], and T = 300 K [blue (dark gray) line]. The temperature 
offset is a = 10%. 

We next examine what happens to the current in the 
presence of a small on-site quadratic potential with on- 
site spring constant ko = 0.001 eV/A 2 u. Shown in Fig. [5] 
are plots of the current as it evolves in time. Compared 
to the current shown in Fig. [31 although there is still 
the initial overshoot to a negative value, the current ap- 



pears to approach a quasi-steady-state value. However, 
the value of the on-site potential is still not enough to 
suppress the large current oscillations, especially when 
the temperature is high. 




150 200 250 
time [10" 14 s] 



FIG. 6. (Color online) Plots of the current when the on-site 
spring constant is ko = 0.1 eV/A 2 u. The left and right chains 
have lengths (a) 7V L = N R = 100 and (b) N L = N R = 50. 
The average temperatures between the chains are T — 100 K 
[green (lowest gray) line], T = 300 K [red (middle gray) line], 
and T — 500 K [blue (upper gray) line]. The temperature 
offset is a — 10%. The dash lines are the values of the steady- 
state current, corresponding to T = 100 K, 300 K, and 500 K, 
calculated independently from the Landauer formula. 

Based on our results for the small on-site ko, we ex- 
pect the current to show a more prominent approach to 
a quasi-steady-state value when we increase the value of 
the on-site ko further. Shown in Fig. [5] are plots of the 
current when ko — 0.1 eV/A 2 u, i.e., at 10% of the value of 
k. Also shown in the figure are dashed lines representing 
the steady-state current calculated using the Landauer 
formula. We now see that the current approaches a quasi- 
steady-state value and that this quasi-steady-state lasts 
longer for longer leads. Note, however, that the quasi- 
steady-state lasts for more than 2t m (t m is defined to 
be the time when sound waves travel the left or right 
chain). After time 2t m , the waves or disturbances that 
have been reflected back at the hard walls at the edges of 
the leads have returned back to the interleads coupling 
and interfere with the other waves there. This results in 
the current beginning to oscillate wildly. 

We next set the leads to have the same temperature 
and, therefore, according to the Landauer formula, we 
should not expect to have current flowing in the long- 
time steady-state limit. We do want to know if we get a 
nonzero transient current and, if so, how would it behave 
even when there is no temperature bias. Shown in Fig. [7] 
are plots of the current flowing out of the left lead when 
the leads have the same temperature and there is no on- 
site potential. Also shown are plots of the total energy 
that has flowed out of the left lead. The total energy at 
time t is calculated by taking the area under the curve 
for the current up to time t. Since the leads are indistin- 
guishable, the current plots for the left and right leads 
are exactly the same. 

Energy is added to the system when the leads are ini- 
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FIG. 7. (Color online) Plots of (a) the current and (b) the 
energy as functions of time when fco = 0. The left and right 
chains have the same temperature, i.e., a = 0, with tempera- 
tures T = 10 K (black line), T = 100 K [red (gray) line], and 
T = 300 K [blue (dark gray) line]. The chains have length 
N L = N R = 50. 



tially attached. This energy flows into both the left and 
the right leads. So even when there is no temperature 
bias, because of the externally added energy, a thermal 
current appears in the system. However, because there is 
no on-site pinning potential, the current does not settle 
into a quasi-steady-state value. The phonons are elasti- 
cally bouncing back and forth between the fixed walls at 
the edges of the two leads. 
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FIG. 8. (Color online) Plots of (a) the current and (b) the 
energy as functions of time when the on-site spring constant 
is fco = 0.001 eV/A 2 u. The left and right chains have the 
same temperature of values T = 10 K (black line), T = 100 K 
[red (gray) line], and T = 300 K [blue (dark gray) line]. The 
chains have length JVl = ATr = 50. 

Shown in Fig. [5] are plots of the current and to- 
tal energy when a small on-site potential with fco = 
0.001 eV/A 2 u is present in the system. The presence 
of the on-site pinning potential suppresses the current to 
a quasi-steady-state value. The total energy shown in 
Fig. ISJb) also displays a tendency to settle down to a 
quasi-steady-state value. 

Shown in Fig. [9] are plots of the current and total 
energy when the on-site potential has spring constant 
fco = 0.1 eV/(A 2 u). The value is large enough for the cur- 
rent to settle into a quasi-steady-state value. However, 
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FIG. 9. (Color online) Plot of the energy as a function of time. 
The inset shows the plot of the current as a function of time. 
The left and right chains have the same temperature T = 
300 K and the on-site spring constant is fco =0.1 eV/A 2 u. 
The chains have length Al = ATr = 50. 



since the leads are finite, the system only has limited time 
before the reflected phonons begin to arrive and interfere 
with the other phonons moving through the interleads 
coupling. 

We can also calculate the current at sites away from 
the midpoint of the composite system where the sudden 
attachments occurred using Eq. (fl~5|) . For a composite 
lead containing Ni, = 100 and 2Vr = 100 sites, we deter- 
mine the current at sites n\ = TVl + 50 and = A^l + 75. 
Shown in Figs. [TU] and [TT] are plots of the current at those 
locations. In Fig. [TUl the spring constant is 1 eV/(A 2 u) 
and there is no on-site potential. Because n\ and n-i are 
away from the midpoint, the current is going to take some 
time, depending on how fast phonons travel in the chain, 
to reach them. Defining a time unit [t] = 10~ 14 s, it takes 
50 [t] for the current to reach n\ and 75 [t] to reach 712, 
therefore implying a speed of 1 site per [t]. 
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FIG. 10. (Color online) Plots of the current at various loca- 
tions on the right lead. The length of the left and right leads 
are jVl = -/Vr = 100. The darker (black) lines in both panels 
(a) and (b) are the current at A^l- The current at m = ^+50 
[blue (dark gray) line in panel (a)] and 712 = A^l + 75 [red 
(gray) line in panel (b)] are also shown. The on-site spring 
constant is fco = and the temperature T = 300 K. 

In Fig. [TT] a quadratic on-site potential with spring 
constant fco = 0.1 eV/(A 2 u) is added to the system. The 
two leads have the same temperature T — 300 K. For the 
midpoint where the sudden connection at t = occurs, 
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upon connection the current immediately drops and then 
oscillates and decays until it reaches a quasi-steady-state 
value which, for this case, is zero. Away from the mid- 
point, the current takes some time to reach the site and 
so initially we see a flat line at zero. Since ko is not zero, 
the speed of the current is not exactly 1 site per [t] as we 
found in Fig. [TUJ 




3I)() 



time [10 



FIG. 11. (Color online) Plots of the current at different loca- 
tions on the right lead when N L = Nn = 100, T = 300 K, and 
on-site spring constant ko = 0.1 eV/A 2 u. Fhe darker (black) 
lines in both panels (a) and (b) are the current at Nl. Fhe 
current at n\ [blue (dark gray) line in panel (a)] and ri2 [red 
(gray) line in panel (b)] are also shown. 



and the initial chain temperatures. A quasi-steady-state 
current is established earlier when the on-site pinning po- 
tential is stronger. Furthermore, when the quasi-stcady- 
state current is established, its value turns out to be the 
same as the steady-state value calculated using the Lan- 
dauer formula. Although in the standard modeling of 
bosonic heat baths infinite systems are usually employed, 
we show here that integrable finite systems also behave 
like an infinite system but only within short time scales 
proportional to the system size N and provided that a 
small on-site pinning potential is present. At time scales 
longer than N, we have recurrent and quasi-periodic be- 
havior. It is interesting to investigate what happens if 
the chains are nonlinear; however, the eigenmode expan- 
sion technique discussed in this paper is unable to handle 
nonlinear systems. 
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Appendix A: Recovering the Landauer formula 



V. SUMMARY AND CONCLUSION 

We examine the energy current in a quantum system 
consisting of two finite chains that are abruptly attached 
at time t = 0. We calculate the current using an eigen- 
mode approach that makes use of a similarity transfor- 
mation extending the normal modes of a single harmonic 
oscillator into a many-particle, but finite, chain of har- 
monic oscillators. In addition, we use a nonequilibrium 
Green's function approach to verify some of the results 
from our normal mode calculations. We find that the 
presence of a quadratic on-site pinning potential is cru- 
cial in establishing a quasi-steady-state current. In a fi- 
nite system where there is no such on-site potential, the 
current does not settle to a quasi-steady-state value, even 
when the two original chains have the same temperature 
at a value as low as 10 K. The phonons would simply 
bounce back and forth elastically between the two fixed 
walls at the edges of the chains. In the presence of the on- 
site pinning potential, we find a tendency for the current 
to establish a quasi-steady state (also, see the Appendix). 
This crucial role of the on-site pinning potential in estab- 
lishing a steady-state current should also be present even 
when the leads arc semi-infinite. Computationally, the 
presence of the quadratic on-site pinning potential ren- 
ders well-behaved Green's functions. The form of the 
time evolution of the current depends on the interplay 
between the strength of the potentials, the chain length, 



In this appendix, we simplify the formula for the cur- 
rent based on Eqs. (|12| to (fl"4|) and derive the Landauer 
formula in the large-size limit. We first introduce a new 
notation and then group the terms into two parts. The 
first part involves terms that are proportional to the dif- 
ference between the Bose distributions of the left and 
right leads. The rest of the terms constitutes the other 
part and are interpreted to be the transient contribution 
to the current. We assume that each chain has length N 
and, therefore, that the total size of the composite chain 
is 2N. Define 



k) K 



N 

E (fc|ni)<«i|*>, 

ni — 1 
2JV 

E (*|ni)<ni-JV|fe> 
m=N+i 



(Al) 



where (n \k ) = ( k\ n) 
■gL, (k\n) = (n\k) = 



Tj- sin(fcn), 



7TJ 

2N+1 ■ 



Notice that (k \k ) R 



2 pf +1 sin(fcn), and k — 
= (k\ k) L (-l) J+] and 



( k\ k 2 ) L (k 2 \N + 1 > = < k\ k 2 ) L (k 2 \N ) The 
expression for the current I L (t) can be separated into 

^stdy ( 



two parts, I^ raas (t) and ii dv (t), where the transient con- 
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tribution is 

i„ s (^H 2 E(^ + 4 R + 1 ) 

k 



x £ cosK 2 i) (fc| fc 2 ) L (fc 2 |iV) 

A; 2 odd 

x E 



0J k ' 



k± even 

t)(fc|fc 2 ) L (fc 2 |AO 

&2 even 

and the steady-state contribution is 



(A2) 



E 

fci even 



x ^ cos( Wfc2 <)(fc|fc 2 ) L (fc 2 |A) 

&2 odd 

^o 2 E(^-/f) 



x E cos(wfe^) (A^|fci)(fci|fc 

fci even 

X 



k2 odd 



E — + — sin( Wfc2 <) (fc|fc 2 ) L (fc 2 |A). 

(A3) 



In the above, the summation for fc extends over all ]^fj 
for j — X,...,N, the summation involving "fci even" 
r ■ } 2N 

is on ti G i 2N+1 f i Ji is even, and the summa- 

tion involving "fc 2 odd" is on fc 2 e j j , j 2 

is odd, and so on. The dispersion relation satisfied is 
uj q = \j2k(\ — cos q) + fco where = fc, Wi = fco, 
/? = l/fe^^s - 1), and a = L, R. The sum in Eq. (|AT]) 
can be carried out analytically, resulting in 

-1 



(iV|fci)(fci|fc) L 



(27V + 1)VW+1) 



x sin(fciV) 



JL 



cos fci — (—1) 
cos k — cos fci 



(A4) 



A similar expression can be derived for ( k\ k 2 ) (k 2 \N). 
Using Eqs. (|A2[) to (| A4|) . the current can now be calcu- 
lated in computer time proportional to 0(N 2 ). This is 



in contrast to NEGF calculations which go as 0(N 3 ) in 
computational complexity. 

All of the expressions derived at this point are exact. 
We now make an approximation in order to extend our 
calculations for the steady-state contribution to large N 
and eventually arrive at the Landauer formula. Notice 
that in Eq. (|A4[) and the corresponding expression for fc 2 
that the terms involving ki rs k « fc 2 would dominate 
the summation, especially when N approaches infinity. 
Consequently, 



J s L tdyW ~M) 



1 



(fl-fl)^k 



1 1 ^ sin (w£ - Ld k2 ) t 



2N +12 t—-' cos k — cos fc 2 

k 2 z 



— (e 

27V + 1 1 ^ 

kfci=e 



-1 



(A5) 



cos k — cos fci 



E 

fci odd 



l 



cos k — cos fci 



Let N — > oo and then followed by f — > oo, we get 

1 1 ^ sin (wg - u k2 ) t ^ 1 
2iV+12^ cosfc-cosfc 2 ~ 2sinfc' 

Furthermore, we have 



lim 

tv->oo 2A 



(A6) 




1 



(A7) 



fci odd 



cos k — cos fci I sin fc 



for some fc € | | . . We then recover the Landauer 
formula 



J stdy 



1 



1 



N ■ 



1 ^ 

k 



Ul: /f)sinfc 



1 

2^ 



(A8) 



Note in the above that the discrete summation over wave 
vector fc is converted into a continuous integration over 
the angular frequency uj. Furthermore, we want to em- 
phasize that although the on-site constant fco appears 
in the expressions for both / s tdy and /trans, its presence 
in the steady-state contribution / s tdy does not prevent 
the contribution to approach a steady-state value in the 
long-time limit. In contrast, the value of fco is crucial for 
the transient contribution /trans to decay away. A zero 
fco would result in /trans having a strong time-dependent 
zigzag-like behavior that would dominate the total energy 
current at all times. There would be no steady current 
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flow even when TV — > oo. However, even a small on-site contribution to decay away in the long-time limit and 
potential, say ko/k — 0.1, would result in the transient leave only the contribution from the steady-state term. 
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